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Abstract 

The  effects  of  prebuckling  displacements  on  the  buckling  of 
laminated  composite  circular  cylindrical  shells  are  investigated. 

Both  axial  compression  and  pure  torsion  are  considered  for  two  shell 
geometries.  A  clamped  prebuckling  boundary  condition  is  used  for  all 
analysis  with  four  buckling  boundary  conditions  applied  during  the 
buckling  process.  The  shell  walls  are  made  up  of  a  6  ply  laminate 
with  several  symmetric  ply  orientations.  The  study  was  made  using 
the  STAGS  computer  code,  utilizing  the  linear  bifurcation  branch  with 
linear  prebuckling  displacements.  The  results  are  compared  to  the 
buckling  loads  determined  when  prebuckling  displacements  are  neglected. 
It  is  shown  that  prebuckling  deformations  generally  tend  to  decrease 
the  buckling  load  of  a  composite  shell.  Increased  buckling  loads  can 
occur  under  axial  compression  with  prebuckling  displacements  assumed 
present,  for  particular  ply  orientations,  due  to  a  higher  bending 
stiffness.  Similarly,  under  torsion,  an  increase  in  buckling  load 
can  occur  because  of  a  higher  tension  hoop  stress.  It  is  also  shown 
that  prebuckling  displacements  can  cause  shell  buckling  before  failure 
of  the  fibers  occurs. 
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BUCKLING  ANALYSIS  OF  LAMINATED  COMPOSITE 


CIRCULAR  CYLINDRICAL  SHELLS 


I.  Introduction 


Previous  Work 

Laminated  composite  materials,  due  to  their  high  strength  and 
lightweight  properties,  are  rapidly  replacing  metals  in  many  structural 
applications.  A  widely  used  structural  element,  the  circular  cylin¬ 
drical  shell,  is  commonly  maoe  of  various  metals  tor  such  applications 
as  missile  interstage  and  intertank  structures;  missile  payload, 
guidance  and  control  adapters;  satellite  components;  and  aircraft. 
Studies  have  shown  Ref  (1)  that  composites  reinforced  by  advanced 
high  modulus  fibers  such  as  boron  and  graphite  permit  weight  savings 
of  25-40%  over  competing  metallic  designs.  These  advanced  composite 
materials  are  being  used  or  considered  for  use  in  many  military  and 
civilian  vehicles  such  as  the  F-lll,  F-14,  F-15,  F-16,  F-18,  C-5A 
and  space  shuttle. 

With  the  high  strength  and  lightweight  properties  of  composites, 
it  can  be  seen  that  the  detail  study  of  composite  cylindrical  shells 
is  important.  One  major  area  of  concern  with  any  structural  element 
Is  with  buckling.  Generally,  composite  shells  will  have  post-buckling 
strength  and  integrity,  however  excessive  buckling  can  lead  to  failure 
due  to  exceeding  maximum  strain  requirements  or,  if  loads  are  repeated, 
exceeding  the  fatigue  life.  There  has  been  a  great  deal  of  work  done 
studying  the  buckling  of  isotropic  cylindrical  shells.  Investigations 
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of  how  the  bucklirg  strength  is  altered  due  to  variations  in  the 


diameter  to  thickness  ratio,  length  to  radius  ratio,  initial  imper¬ 
fections,  prebuckling  deformations  and  buckling  boundary  conditions 
for  a  variety  of  loading  conditions,  has  been  done.  For  a  complete 
list  of  references  on  buckling  of  isotropic  shells,  see  Reference  (2). 

Jones  and  Hennemann  (1)  invescigated  the  effects  of  prebuckling 
deformations  on  buckling  of  laminated  composite  circular  cylindrical 
shells.  For  a  variety  of  antisymmetric  cross  ply  orientations,  they 
looked  at  a  simply  supported  edge  boundary  condition.  The  effects  of 
prebuckling  displacements  were  found  to  be  negligible  for  buckling 
analysis,  and  they  suspected  similar  results  for  all  other  boundary 
conditions.  Further  work  done  on  stability  of  composite  shells  can 
be  found  in  References  (3)-(9). 

Problem  Definition 

The  purpose  of  this  thesis  is  to  investigate  how  prebuckling 
displacements  effect  the  buckling  of  composite  shells  for  two  length 
to  radius  ratios.  Each  geometry  is  studied  first  under  axial 
compression  and  then  torsion.  Along  with  the  geometric  instability, 
a  second  area  investigated  is  material  instability.  Different  ply 
orientations  for  6  ply  laminates  are  studied  to  compare  buckling  loads 
to  composite  strengths  as  estimated  from  maximum  strain  criteria.  The 
shell  dimensions  considered  in  this  thesis  are  typical  of  current 
laboratory  tubular  specimens  presently  being  tested. 

Using  classical  Fourier  analysis  in  conjunction  with  Flugge's 
shell  equations  modified  for  anisotropic  laminated  materials,  Whitney 
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and  Sun  (10)  have  obtained  buckling  loads  for  the  same  shells 
considered  in  this  thesis.  The  development  of  these  equations  is 
shown  in  Appendix  A.  Their  results  suggest  that  buckling  is  a 
potential  problem  for  tubular  specimens.  This  thesis  takes  into 
account  prebuckling  deformations  that  Whitney  and  Sun  neglected. 

Scope 

This  thesis  looks  at  two  composite  circular  cylindrical  shells, 
under  two  independently  applied  loads,  axial  compression  and  pure 
torsion.  Both  shells  have  a  6  ply  laminate  with  (+  0,  0°)s  orientation, 
illustrated  in  Figure  1.  Computer  runs  are  made  for  6  ranging  from  0° 
to  90°.  The  buckling  boundary  conditions  are  clamped  and  summarized 
as  follows: 


u 

V 

w 

X 

CC1 

FREE 

FREE 

0 

0 

CC2 

0 

FREE 

0 

0 

CC3 

FREE 

0 

0 

0 

CC4 

0 

0 

0 

0 

where  u,  v  and  w  are  displacements  in  the  axial,  circumferential  and 
radial  directions  respectively,  as  shown  in  Figure  2.  The  prebuckling 
boundary  condition  used  for  all  cases  is  CC1.  The  linear  bifurcation 
branch  of  the  STAGS-C  computer  code  is  used,  incorporating  linear 
prebuckling  deformations. 
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II.  Theory 


Bifurcation  Buckling 

Since  the  linear  bifurcation  branch  of  STAGS  is  used  in  this 
investigation,  it  is  appropriate  to  first  discuss  the  concept  of 
bifurcation  buckling.  In  general,  buckling  can  be  defined  as  a 
structural  deformation  due  initially  to  instability  under  load 
(Ref  11) .  It  does  not  matter  whether  the  deformation  is  elastic 
or  permanent,  or  whether  it  leads  immediately  to  collapse  or  not. 

For  a  bifurcation  type  of  buckling,  a  structure  is  said  to  rapidly 
move  from  one  equilibrium  state  to  another. 

Consider  the  case  of  the  cylindrical  shell,  subjected  to  a 

load  P.  The  load  increases  monotonically  from  an  initial  value 

of  zero.  Assume  that  for  sufficiently  small  values  of  P,  there 

exists  a  unique  solution  to  the  shell's  equilibrium  problem.  The 

equilibrium  configuration  given  by  this  solution  is  called  the 

primary  state  or  prebuckled  state.  Then  take  some  value  of  P,  say 

P  *  P^,  where  there  exists  another  equilibrium  configuration  that 

is  infinitesimally  close  to  the  primary  state.  The  existence  of 

two  adjacent  infinitesimally  close  equilibrium  configurations  at 

the  same  value  of  the  load  P  =  F  is  called  bifurcation,  and  P 

b  b 

is  called  the  bifurcation  buckling  load.  This  can  be  illustrated  by 
the  curve  in  Figure  3,  in  which  P  is  plotted  against  lateral  displace¬ 
ment  <5.  At  each  point  on  the  curve  an  equilibrium  configuration  exists. 
The  bifurcation  point  is  the  point  of  intersection  of  the  primary  path, 
which  starts  at  the  origin  and  is  initially  stable,  and  the  secondary 
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LOAD  P 


Figure  3.  Load-Displacement  Curve  for  Bifurcation 


path  at  P  “  P^.  Since  the  stability  on  the  primary  path  is  lost 
above  the  bifurcation  point,  the  shell  will  either  deform  into  the 
equilibrium  configuration  on  the  secondary  path,  or  it  will  be  set 
into  motion.  The  new  equilibrium  configuration  on  the  secondary 
(post  buckled)  path  is  called  the  buckled  state  Ref  (12). 

The  bifurcation  buckling  load  P^  can  be  determined  by  different 
methods.  STAGS  uses  finite-difference  approximations  of  energy 
equations  and  solves  an  eigenvalue  problem.  This  procedure  will  be 
described  in  the  section  that  follows. 

Theory  of  STAGS 

STAGS  (Structural  Analysis  of  General  Shells)  is  a  computer  code, 
developed  to  analyze  general  shells  under  various  static,  thermal  and 
mechanical  loading.  The  approach  used  by  STAGS  is  a  two-dimensional 
finite  difference  approximation  of  energy  relations.  For  this  thesis, 
the  problem  is  for  a  circular  cylindrical  shell,  with  laminated 
composite  fibers  oriented  in  a  number  of  directions.  The  applied 
loads  are  static,  and  thus  the  case  of  static  equilibrium  must  be 
solved.  A  necessary  condition  for  static  equilibrium  is  that  the 
total  potential  energy  V,  must  be  stationary,  which  means  the  first 
variation  of  V  is  equal  to  zero.  The  total  potential  energy  of  the 
shell  is  equal  to  the  strain  energy  U,  of  the  shell  minus  the  work 
done  by  the  external  forces  W, 

V  -  U  -  W  (1) 

where 

W  -  X**  F 
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X*  is  the  vector  of  displacement  components  and  F  is  the  vector  of 
external  forces. 


The  strain  energy  for  an  orthotropic  cylindrical  thin  shell  is 
given  by 


u-i 


Ks  +  a9ce  +  tx6yX0>  dxdy 


(3) 


where  the  strain  terns  ex,  e0  and  yx  are  given  by 

A0 


e„  ■  u,„  +  —  a  2 


‘»x 


2  $x^ 


v’0  +  W  X  1  r  ’ 
e0  *  - a -  +  7  6e 


&x  =  -w,x 

W,0 

Se  “  -  a 


(4) 


U, 


Txe 


<•  +  »•*>  + 


The  quantities  u,  v  and  w  are  displacement  components  in  the  x,  0  and 
directions  respectively,  and  "a"  is  the  undeformed  middle  surface 
radius.  The  stress  terms  are  found  by 

[>j  ]  *  [Qy  1  [ej  ]  i.J  -  1.2,6  (5) 

where  e. ,  e~,  e,,  on ..  c,,,  and  a,  are  defined  as  e  ,  e  ,  y  ,  a  ,  a 
1  2  6’  1  2’  6  x  0  x0  x  0 

and  xX0  respectively.  The  matrix  is  called  the  transformed 
reduced  stiffness  matrix  (see  Ref  (13)  for  more  detailed  analysis 
of  Qy  )  . 

For  the  case  of  the  angle  ply  laminate,  the  stiffness  matrix  will 
be  different  for  each  layer.  It  is  therefore  necessary  to  look  at 
the  resultant  forces  and  moments  acting  on  the  laminate.  These  are 
obtained  by  integrating  the  stresses  in  each  layer  through  the 
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N-layered  laminate  thickness  as  follows 

.t/2  N 

N  *  o.dz  =  l  Ojdz  1  =  1,2,6 

-t/2  1  K=1  JZR_L 

and 

(t/2  N  fZK 

M.  =  Za.dz  =  l  Za.dz  i  =  1,2,6 

-t/2  K=1  JzK_i 

Using  the  relations  for  the  change  in  curvatures: 

ge,e 

Kx  *  ^x,x  Kd  a 


< 


x9 


°x ,  6 


+  6 


9,x 


) 


(6) 


(7) 


(8) 


and  combining  this  with  Eqs  (4)  and  (5) ,  Eqs  (6)  and  (7)  become  in 
matrix  form 


i,j  -  1,2,6 


where  the  stretching,  bending-stretching  coupling,  and  bending 
stiffness  submatrices  are  defined  as 


(Ay,  B^j ,  D±j)  = 


t/2 


J-t/2 


(1,  Z,  Z2)  Qjjdz 


(9) 


(10) 


For  the  shell  that  is  laminated  symmetrically  with  respect  to  the 
middle  surface,  as  is  the  case  in  this  thesis,  the  bending-stretching 
submatrix  Bjj  equals  zero.  This  yields  the  following  expressions  for 
middle  surface  strains 

<ej°> a  CAij  l"1  <Ni>  <n) 
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and  curvatures 


{<d  >  =  n  D±J  □~1  {M±>  (12) 

With  the  middle  surface  strains  and  curvatures,  the  total  strains  at 
an  arbitrary  distance  z  from  the  middle  surface  are 

GJ  "  £j°  +  ZKj  i  *  1»2’6  (13) 

Using  Eqs  (13)  and  (5) ,  all  the  terms  in  the  strain  energy  are  now 
known.  In  order  to  approximate  the  strain  energy  and  thus  the  total 
potential  energy,  STAGS  models  the  shell  structure  by  placing  mesh 
lines  paralled  to  the  x  and  y  coordinate  lines.  The  intersection  of 
these  mesh  lines  represent  nodes  or  mesh  stations  i,  and  the  area 
bounded  by  mesh  lines  would  be  a;  the  area  of  the  i1"*1  subregion.  The 
total  strain  energy  would  then  be  given  by 

M  ,  . 

U  *  £  A  u1  ’  a1  (14) 

i 

where  M  is  the  number  of  mesh  stations  and  Au*,  which  is  the  integrand 
of  Eq  (3),  is  the  strain  energy  density  at  mesh  station  i.  With  the 
work  W,  and  strain  energy  U,  now  given  in  terms  of  the  unknown 
displacements,  STAGS  replaces  finite  difference  approximations  for 
those  displacements  in  the  total  potential  energy  V,  given  by  Eq  (1). 

As  previously  stated,  for  the  condition  of  stationary  potential 
energy,  the  first  variation  must  vanish.  This  leads  to  the  equation 

LX  -  F  (15) 
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where  L  is  a  stiffness  operator  which  relates  displacement  components 
and  external  forces,  and  is  defined  as 

LX  =  GRAD  U  (16) 

Since  only  the  linear  terms  are  included  in  the  deformations  of  strain 
and  curvature,  then  for  linear  bifurcation  theory,  L  is  a  linear 
operator  represented  in  matrix  form. 

For  the  bifurcation  problem  let  Xq  be  a  solution  of  Eq  (15)  for 
a  given  external  force  vector  F.  If  another  displacement  vector  Y, 
in  every  neighborhood  of  Xq  satisfies  the  equation 

1,Y  =  F  (17) 

then  bifurcation  is  said  to  take  place  for  the  shell  under  load  F. 

Since  Eq  (15)  can  have  multiple  solutions,  a  necessary  condition  for 

bifurcation  required  that  L'Y  ,  which  is  the  derivative  of  the  stiffness 

Ao 

operator,  be  a  singular  matrix.  This  leads  to 

det  (L’  )  =  0  (18) 

Ao 

from  which  classical  bifurcation  theory  may  be  obtained.  It  is  then 
assumed  that  XQ  may  be  written  as 

XQ  -  A  XL  (19) 

where  X^  is  the  linear  solution  for  a  load  vector  F^.  Thus  Eq  (18) 
becomes 


det  (L' 


0 


(20) 
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Eq  (20)  represents  an  algebraic  eigenvalue  problem.  Omitting  the 
nonlinear  prebuckling  rotations  for  linear  bifurcation  theory, 

Eq  (20)  becomes 

det  (A  -  XD)  =  0  (21) 

this  leads  to  the  following  eigenvalue  problem 

AX  =  XDX  (22) 


where  A  is  the  linear  stiffness  matrix,  and  D  is  the  bending  stiffness 

matrix.  Briefly,  the  formation  of  the  A  and  D  matrix  is  as  follows. 

Consider  the  stiffness  operator  derivative  L'Y  ,  which  represents  the 

Ao 

coefficient  matrix  for  the  linear  system  of  governing  algebraic 

equations  where  the  nodal  displacements  are  the  unknowns.  The  elements 

of  L'v  are  given  by 
Ao 


L' 


a2u 


(23) 


Applying  Eq  (23)  to  Eq  (14)  yields 


32U  =  V  ak  92Aljk 

3Xi3Xj  k=l  3X^3Xj 

The  strain  energy  density  AU  can  be  written 

AUk  -  ZkTDkUk 


(24) 


(25) 


where  the  matrix  Z  represents  the  strains  and  curvatures  at  mesh 
station  k  and  D  is  a  matrix  of  constants  that  depend  on  the  material 
properties.  Recall  the  integrand  of  Eq  (3)  is  the  strain  energy 
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density  for  the  cylindrical  shell.  The  vector  of  linear  stress 
resultants  at  station  k  is  defined  as 


(26) 


Incorporating  Eqs  (25)  and  (26)  together,  the  k 


th 


term  of  Eq  (24)  is 


32AUk 
3  3 

X,*  X 


27^ 


3^Z 


xi  xj 


kT  V 

+  9Z  Dk  3ZX 

3y 

Xi  Xj 


(27) 


For  the  case  of  neglected  prebuckling  notations,  the  last  term 

of  Eq  (27)  generates  the  linear  stiffness  matrix  A.  The  first  term  of 

Eq  (27)  contributes  to  the  bending  stiffness  matrix  D.  The  critical 

load  F  is  determined  when  the  smallest  eigenvalue  X  is  found  that 
CK. 

satisfies  Eq  (21).  The  critical  load  is  a  multiple  of  the  linear 
force  vector  FL 

rGR  '  >FL  <28) 

Ref  (14)  shows  a  more  detailed  analysis  of  the  eigenvalue  problem. 
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III.  Modeling 


Boundary  Conditions 

There  are  four  clamped  boundary  conditions  used  in  this  thesis 
referred  to  as  CC1,  CC2,  CC3  and  CCA  (see  Fig  A  for  a  physical 
representation).  Rollers  allow  freedom  of  movement  while  flanges 
restrict  displacement  perpendicular  to  the  plane  of  the  flange. 

Due  to  the  restraining  moments  of  a  clamped  condition  the  slope 
w,x  is  zero  at  the  shell  boundary.  Notice  also  that  the  displacement 
w,  is  zero  at  the  boundary  for  all  four  cases. 

STAGS  models  these  four  boundary  conditions  by  allowing  or 
restraining  the  particular  displacements.  The  Nx  =  0  and  Nxy  =  0 
conditions  are  specified  by  allowing  u  and  v  respectively  to  be  free. 
The  u  =  0,  v  =  0,  w  =  0,  and  w,x  =  0  restraints  are  specified  directly 
by  setting  the  appropriate  edge  quantity  to  zero.  STAGS  also  has 
several  built-in  boundary  conditions.  One  of  them  was  used  in  this 
thesis,  referred  to  as  symmetry.  Symmetry  means  that  a  boundary  lying 
along  the  y  or  circumferential  direction  has  u  =  w,x  =  0,  and  a 
boundary  lying  along  the  x  or  axial  direction  has  v  =  w,y  =  0. 

See  Fig  5. 

Basic  Shell  Models 

As  previously  stated,  there  are  two  L/R  ratios  studied  in  this 
thesis.  Each  shell  geometry  has  a  radius  (R)  of  5.25"  and  a  thickness 
(h)  of  .5".  The  L/R  ratios  considered  are  12  and  7,  yielding  a  length 
(L)  of  63"  and  36.75"  respectively.  For  each  shell  there  are  two 
loading  arrangements  applied  independently,  axial  compression  and 


w  =  0  w,x  =  0  Nx  =  0  v  =  0 


torsion.  Due  to  the  nature  of  the  loads,  the  shell  model  for  each 
load  was  different  and  will  now  be  presented. 

For  the  case  of  axial  compression,  a  constant  line  load  was 
applied  to  the  shell.  Due  to  the  symmetry  of  the  loaded  shell,  only 
1/8  of  the  shell  was  modeled.  The  only  difference  between  the  model 
for  L/R  =  12  and  L/R  =  7  is  the  value  of  L/2  (see  Fig  5).  The 
boundary  condition  for  each  boundary,  as  seen  in  Fig  5,  is  as  follows: 
For  boundary  1,  which  is  the  top  of  the  shell,  the  conditions  are 
specified  according  to  the  particular  boundary  condition  being 
investigated.  All  prebuckling  boundary  conditions  are  CC1,  and 
buckling  boundary  conditions  are  CC1,  CC2,  CC3  and  CCA.  Boundaries 
2,  3,  and  A  are  considered  to  be  symmetry. 

For  the  case  of  torsion,  due  to  the  fact  that  symmetry  of 
displacement  is  not  present,  the  whole  shell  had  to  be  modeled 
(see  Fig  6).  A  problem  arose  that  indicated  a  rigid  body  mode,  and 
thus  it  became  necessary  to  improve  the  model.  It  was  determined 
that  fixing  the  circumferential  displacement  v,  of  at  least  one  node 
on  a  boundary  would  alleviate  the  rigid  body  problem.  Some  sample 
runs  indicated  that  holding  all  the  nodes  of  one  edge  fixed  in  the 
v  direction  made  for  a  better  model  than  just  fixing  one  node.  This 
was  taken  care  of  for  that  edge  by  letting  v  =  0  in  the  boundary 
conditions.  Also,  no  shear  for  that  edge  was  applied  since  it  was 
fixed  in  the  circumferential  direction,  which  yields  the  same  as 
applying  the  shear  (see  Appendix  B) . 

For  L/R  =  7,  the  shear  was  applied  to  the  top  edge.  The  boundary 
conditions  for  the  top  and  bottom  edge  were  the  same  conditions  used 
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in  axial  compression,  CC1  prebuckling  and  CC1  through  CC4  buckling. 

The  only  difference  being  the  added  constraint  of  v  =  0  along  the 
bottom  edge.  What  were  called  boundaries  2  and  4  for  the  1/8  shell 
model,  are  not  really  boundaries  for  the  whole  shell.  There  is  an 
option  in  STAGS  that  allows  two  boundaries  to  have  compatible 
displacements,  since  there  must  be  an  allowance  for  four  boundaries. 

This  option  takes  care  of  the  treatment  of  the  whole  shell. 

The  L/R  =  12  geometry  had  the  same  model  with  the  exception  of  a 
different  mesh  size.  The  mesh  sizes  will  be  discussed  in  the  section 
that  follows. 

Both  shells  considered  were  made  of  graphite/epoxy  composites, 
with  the  following  material  properties 

El  »  20  x  106psi,  Ex  =  106psi,  G^T  =  0.6  x  106psi,  vLT  -  0.25 

where  L  denotes  the  direction  parallel  to  the  fibers,  T  denotes  the 
direction  transverse  to  the  fibers,  and  is  the  Poisson  ratio  as 
determined  from  a  tensile  test  parallel  to  the  fibers  (Ref  10). 

The  laminates  for  each  shell  were  6  ply  and  oriented  by  (+0,  0°)s. 
0  was  0°,  15°,  30°,  45°,  60°,  75°  and  90°. 

Finite  Difference  Mesh  Arrangements 

There  were  four  different  mesh  arrangements  used  in  this  thesis. 

For  the  1/8  shell  model,  210  nodes  at  a  mesh  size  of  .92"  x  .92"  were 
used  for  the  L/R  =  7,  and  340  nodes  at  a  mesh  size  of  .92"  x  .95"  for 
L/R  =12.  An  attempt  to  test  the  accuracy  obtained  by  these  mesh  sizes 
was  made  by  tripling  the  number  of  nodes  (which  almost  tripled  the 
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computer  cost).  This  smaller  mesh  size  gave  a  higher  buckling  load 
of  less  than  only  1.5%  which,  for  the  additional  cost,  is  an 
insignificant  increase  in  accuracy.  Decreasing  the  number  of  nodes 
by  30%  gave  only  2  or  3  nodes  per  1/2  sine  wave.  It  was  felt  that 
5  or  6  nodes  per  1/2  sine  wave  were  necessary  for  any  accuracy 
Ref  (15),  and  this  was  achieved  by  the  aforementioned  mesh  arrangements. 

For  the  whole  shell  model,  a  problem  of  ill-conditioned  variables 
arose  for  both  geometries.  The  L/R  =  7  case  was  modeled  with  1054 
nodes  giving  a  mesh  size  of  1.11"  x  1.06".  This  mesh  arrangement  was 
used  for  all  angles  except  6=0°.  At  9  =  0°,  ill-conditioned 
variables  in  the  stiffness  matrix  caused  the  matrix  decomposition  to 
be  discontinued  prematurely  (see  Appendix  B) .  This  was  a  problem 
for  angles  0,  less  than  8°.  At  8°,  the  analysis  was  successful,  and 
from  this  and  other  results,  the  buckling  load  at  9  =  0°  was 
extrapolated.  The  L/R  =  12  model  initially  had  989  nodes  at  a 
mesh  size  of  1.47"  x  1.50",  but  this  mesh  arrangement  yielded  ill- 
conditioned  stiffness  matrix  for  most  all  angles  of  9.  The  number 
of  nodes  was  increased  to  1128  giving  a  mesh  size  of  1.37"  x  1.37". 

The  latter  arrangement  worked  for  all  angles  9  except  0°.  The 
identical  procedure  used  in  obtaining  the  buckling  load  at  9  =  0° 
for  L/R  =  7,  was  used  here  for  L/R  =  12. 
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IV.  Results 


Shell  Buckling  Without  Prebuckling 

The  problem  of  determining  the  buckling  load  for  a  composite 
cylindrical  shell,  without  taking  into  account  the  prebuckling 
deformations,  can  be  done  by  using  a  series  solution.  This  approach 
was  taken  for  an  L/R  =12,  by  Whitney  and  Sun  (10),  for  both  axial 
compression  and  torsion,  and  is  given  in  Appendix  A. 

Whitney  and  Sun  took  the  buckling  loads,  obtained  by  the  series 
solution,  and  normalized  them  for  specific,  ratios  of  L/R  and  R/h. 

This  was  done  by  dividing  the  buckling  load  by  the  product  of 
transverse  modulus  of  elasticity  of  the  fibers  E^,  and  shell  wall 
thickness  h.  This  normalized  load  would  then  hold  for  any  shell  of 
length  L,  radius  R,  and  thickness  h,  provided  the  ratios  of  L/R  and 
R/h  remained  the  same.  The  normalized  loads  were  then  plotted  against 
the  ply  orientations  0.  Included  in  the  graphs  were  an  estimated 
strength  based  on  maximum  strain  criterion,  as  presented  by  Petit 
and  Waddoups  (16) .  Ply  strain  allowables  chosen  for  the  maximum 
strain  criterion  were 

e^+  =  0.01  in/in  eL  =  0.008  in/ in 

e,j+  =  0.004  in/ in  =0.01  in/ in 

Ylt  =  0.015  in/in 

where  +  denotes  tension  and  -  denotes  compression.  The  graphs 
described  above,  for  the  ply  orientation  (+  0,  0°)s,  have  been 
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reproduced  in  Figs  7-9.  Simple  comparisons  can  then  be  made  to  the 
buckling  of  shells  when  prebuckling  deformations  are  included.  The 
following  sections  will  describe  the  effects  of  prebuckling. 

Prebuckling  Effects  For  Compression 

As  previously  stated,  CC1  prebuckling  boundary  conditions  were 
used  in  all  test  cases.  Starting  with  the  shell  having  L/R  =  12, 
runs  were  made  at  each  angle  0  for  CC1  buckling  boundary  conditions. 
The  results  showed  that  for  9=0°,  there  was  a  22%  drop  in  the 
buckling  load  from  the  series  solution  with  no  prebuckling.  For  all 
other  ply  orientations,  the  buckling  load  ranged  from  2%  to  12%  of 
a  decrease  from  the  series  solution.  The  buckling  leads  for  CC3 
buckling  turned  out  to  be  almost  identical  to  the  loads  found  for 
CC1,  all  within  3%.  This  similarity  is  not  unexpected  since  in  both 
cases,  the  axial  displacement  is  allowed  to  move  freely.  Fig  7  shows 
the  comparison  of  the  buckling  loads  determined  from  the  series 
solution,  and  the  buckling  loads  calculated  where  prebuckling  is 
considered.  Appendix  C  shows  the  results  in  tabular  form. 

For  both  CC2  and  CC4  conditions,  the  buckling  loads  were  also 

nearly  identical,  and  for  ply  orientations  except  9  =  15°,  the  loads 

o 

were  the  same  as  for  CC1  and  CC3.  For  0  =  15  ,  however,  an  increase 
in  the  buckling  load  of  approximately  10%  over  the  series  solution  was 
found.  That  is  nearly  20%  higher  than  the  value  expected  from  other 
observed  cases.  Attempts  to  insure  that  the  increase  was  not  caused 
by  a  poor  model  were  made.  Steps  to  check  the  model  included  the 
following;  increasing  the  load  by  50%,  increasing  the  number  of  nodes 
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Figure  7.  Prebuckling  Effect  -  Axial  Compression 
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by  30%  and  modeling  the  entire  shell  with  989  nodes  and  a  mesh  size 
of  1.47"  x  1.50".  In  all  three  cases  the  buckling  load  remained  the 
same  for  9  =  15°. 

Once  it  became  apparent  that  the  model  was  good,  it  became 
necessary  to  explain  why  the  shell,  at  the  angle  6  =  15°,  was  stiff er 
than  for  the  case  of  neglected  prebuckling  displacements.  Investigation 
of  the  prebuckling  moments,  showed  that  for  6  =  15°,  the  values  of 
were  higher  than  for  other  ply  orientations.  With  the  case  where 
no  prebuckling  deformation  is  considered,  the  only  load  is  the 
applied  load  Nx,  and  so  other  prebuckling  forces  and  moments,  including 
are  zero.  Taking  into  account  the  prebuckling  forces  and  moments, 
which  arise  when  prebuckling  displacements  are  considered,  the 
substantially  larger  Mx  for  9  =  15°  makes  for  a  relatively  higher 
bending  stiffness  than  for  other  ply  orientations.  This  relatively 
higher  bending  stiffness,  which  will  resist  buckling  when  the  load 
was  in  the  axial  direction,  makes  for  a  higher  buckling  load  relative 
to  the  case  where  no  prebuckling  displacement  is  considered. 

Test  runs  made  for  the  L/R  =  7  case  were  made  for  the  same 
boundary  conditions  and  ply  orientations.  It  was  found  that  all 
comparisons  made  for  L/R  =  12,  were  similar  to  those  for  L/R  =  7. 

The  only  difference  was  that  the  increased  buckling  load  found  at 
0  *  15°  for  L/R  =  12,  was  true  for  both  0  =  15°  and  0  =  30°  for 
L/R  =  7.  The  prebuckling  ^  was  higher  for  both  angles  as  well,  and 
therefore  the  same  conclusions  are  drawn.  It  was  also  shown  that  the 
change  in  buckling  load  due  to  a  different  L/R  ratio  was  very  small 
(see  Fig  7),  except  at  9  =  15°  and  0  =  30°.  It  can  therefore  be 
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concluded  that  the  boundary  effect  is  very  slight  except  at  6  =  15° 
and  30°. 

It  became  apparent  that  for  both  shell  geometries,  the  significant 
buckling  constraint  is  u  =  0.  It  is  when  the  axial  displacement  on 
the  boundary  is  held  fixed  at  buckling  after  it  is  allowed  to  move 
freely  in  prebuckling  that  the  variation  from  series  solutions  is 
obtained. 

Comparisons  can  also  be  made  between  the  buckling  loads  and 
composite  strength,  as  estimated  from  maximum  strain  criterion.  For 
certain  small  angles  0,  the  curve  for  CC1  and  CC3  buckling  remain  under 
the  maximum  strain  curve,  where  for  the  same  0  for  the  series  solution 
is  above  (see  Fig  7). 

Prebuckling  Effects  For  Torsion 

The  prebuckling  and  buckling  boundary  conditions  used  for  torsion 
were  the  same  as  with  axial  compression.  It  was  found,  as  was  the 
case  with  axial  compression,  that  for  CC1  and  CC3  conditions,  the 
buckling  loads  were  the  same.  Similarly  with  CC2  and  CC4. 

For  cases  of  L/R  =  7  under  torsion  the  values  of  buckling  loads 
shown  in  Fig  8  are  higher  than  the  series  solution,  rang.'ng  from 
15%-28%  (recall  that  the  series  solution  is  done  for  L/R  =  12),  which 
only  indicates  that  a  shorter  shell  demands  a  greater  applied  twist  to 
buckle.  The  shell  length  played  more  of  a  role  in  the  case  of  torsion 
than  it  did  with  axial  compression,  where  very  little  change  in  buckling 
load  occurred  due  to  the  change  in  length  (see  Figs  7  and  3). 

The  final  investigation  was  for  an  I./R  =  12  under  torsion.  It 
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Figure  8.  Prebuckling  Effect  -  Pure  Torsion 
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was  found  that  for  angles  of  6  =  75°  and  below,  the  buckling  loads 
ranged  from  6%  to  12%  below  the  series  solution,  for  CC2  and  CC4 
buckling.  For  CC1  and  CC3  buckling  the  drop  ranged  from  9%  to  22% 

(see  Fig  9).  At  0  =  90°  for  all  four  boundary  conditions,  the 
buckling  load  was  approximately  8%  higher  than  the  series  solution. 
Investigation  of  the  prebuckling  stresses  showed  that  there  is  more 
tension  in  the  hoop  (v)  direction  at  9  =  90°,  than  for  all  other  angles 
of  0.  The  tension  stress,  along  with  the  greater  hoop  stiffness  comes 
about  because  the  torsion  stress  field  acts  over  a  larger  area  in  v 
direction  for  ply  orientation  0  =  90°.  This  is  true  because  the  ply 
orientation  of  90°  is  parallel  to  the  v  direction.  The  tension  stress 
acts  to  resist  torsion  buckling  and  so  this  added  resistance  that  is 
not  as  great  as  for  other  ply  angles  0,  has  the  effect  of  increasing 
the  buckling  load  relative  to  the  series  solution. 

Comparing  the  two  L/R  ratios  for  the  case  of  torsion,  it  is  shown 
that  for  0  >  45°  there  is  a  marked  difference  in  the  shape  of  the 
curves.  For  L/R  =  7,  the  maximum  buckling  load  occurs  at  around 
0  =  60°,  the  curve  then  begins  to  drop  off  as  the  boundary  conditions 
come  in  to  play.  For  L/R  =  12,  the  curve  flattens  out  at  9  >  75° 
because  of  the  increase  in  tension  hoop  stress. 
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V.  Conclusions 


From  the  results  generated  in  this  thesis,  several  conclusions 
can  be  drawn  concerning  the  prebuckling  effects  of  composite  shells, 
with  ply  orientation  of  (+9,  0°)s,  under  axial  compression  and 
torsion. 

1)  Under  axial  compression  and  torsion,  prebuckling  deformations 
generally  tend  to  reduce  the  buckling  load  of  composite  shells. 

2)  For  certain  ply  orientations,  the  prebuckling  effects  can 
increase  the  buckling  load. 

3)  For  a  shell  under  axial  compression,  the  u  =  0  constraint 
in  buckling  is  the  important  factor  in  determining  the  prebuckling 
effects. 

4)  A  change  in  length  of  a  composite  shell  affects  the 
buckling  load  much  more  under  torsion  than  for  compression,  as  is 
consistent  with  isotropic  shells. 

5)  For  both  axial  compression  and  torsion,  the  buckling 
boundary  conditions  CC1  and  C.C3  act  nearly  the  same  while  CC2  and 
CC4  act  similarly. 

6)  For  small  ply  angles  8,  under  axial  compression,  the  effects 
of  prebuckling  displacement  can  cause  buckling  before  failure  of  the 
fibers  occurs,  based  on  maximum  strain  criterion.  This  is  also  true 
under  torsion,  for  angles  0  from  about  25°  to  35°. 

7)  A  higher  prebuckling  bending  stiffness  for  a  particular 
ply  orientation,  under  axial  compression,  has  che  effect  of  increasing 
the  buckling  load  over  that  of  the  solution  generated  with  prebuckling 


8)  A  greater  prebuckling  tension  stress  can  cause,  for  90° 
oriented  fibers,  a  higher  buckling  load  under  applied  torsion,  than 
the  load  determined  when  no  prebuckling  is  considered. 

9)  For  the  case  of  torsion,  the  boundary  conditions  have 
greater  affect  on  the  shell  with  L/R  =  7  than  with  L/R  =  12. 
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Appendix  A 


Series  Solution  Without  Prebuckling 

The  buckling  problem  of  a  composite  cylindrical  shell  was  solved  by 
Whitney  and  Sun  (10).  They  used  classical  Fourier  analysis  in  conjunc¬ 
tion  with  Fliigge's  shell  equations  modified  for  anisotropic  laminated 
materials.  Parts  of  the  development  will  be  shown  here  to  briefly  out¬ 
line  the  procedure.  For  the  complete  analysis,  see  Reference  (10). 

Consider  a  laminated  cylindrical  shell  of  wall  thickness  h  with  a 
middle  surface  radius  of  R.  The  coordinate  system  is  in  the  middle  sur¬ 
face  with  x,  0  and  z  measured  along  the  longitudinal,  circumferential, 
and  radial  direction,  respectively.  Denoting  partial  differentiation 
by  a  comma,  the  governing  equations  for  the  static  stability  analysis  of 
laminated  anisotropic  cylindrical  shells  based  on  Fliigge's  theory  (17) 
are  of  the  form  (18) 

Lnu"  +  hi'’’  +  Ln”  +  N*lu’xx 

N  *  u°  N* 

0  “’00  .  .  -  x9  o 

+  Wx)  +  2  ~  u,x9  =  0 


L12U°  +  L22v°  +  L23w  +  Nx  V’xx 

N01  Nx0 

+  -T  <w’6  +  v’90>  +  2  ~T  (w’x  +  v’xe>  =  0 

K  1 

L, _u°  +  L0„v°  +  L,,w  -  N  w,  -  —  (u’ 
13  23  33  x  ’xx  R  ’x 

v’«  .  w,ee,  „  s*8  .  .  , 

--TT  +  — )  -  2  -jf-  (».x6  -  »,„)  -  o 


where  u°  and  v°  are  the  inplane  displacements  of  the  middle  surface  in 


the  x  and  0  directions,  respectively,  w  is  the  transverse  displacement 

of  the  middle  surface,  and  N  i,  N*  and  N*  are  initial  (pre-buckling) 

x  y  xo 
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membrane  force  resultants  defined  in  the  usual  manner.  Looking  at  three 
of  the  six  operators  ,  gives  an  indication  of  the  complexity  of  the 
problem.  The  L^.  operators  are  defined  for  symmetric  laminates. 

hl  =  All<  >»xx  +  2^<  >’x0+t!<  >*00  (M) 


L13  =  -  -f  (  }’*xx 


h*  (  )  +^6  (  ) 

1?  '  *xx0  r3  ( 


,  °26  ,  .  A12  .  .  1  ,  °26  .  (A5) 
+  <  >  »000  +  “  (  )  >x  +  ^2  (A26  +  (  )  ’0 

L33  =  °11(  }’xxxx  +  4  (  ^  *xxx0  +  ^2  (D12  +  2D66) (  ^’xxee 

+  4%().xeee+%(>.eeeo-2^l(>.xo  <A« 

®22  1  ^22 
+  2  ,(  }  ’60  +  R1  (A22  +~^H  } 

where  A„  and  are  the  inplane  stiffness  matrix  and  bending  stiffness 
matrix,  respectively  which  are  defined  in  Eq. (10).  In  the  present  anal¬ 
ysis,  Initial  loads  are  constant  and  of  the  form  N  ^  =  k,  N  ,  N  *  =  k~N  , 

x  o  0  2  o 

and  =  k_N  ,  where  k.  are  prescribed  constants,  and  N  is  a  charac- 
x:  3  o  l  r  ’  o 

teristic  loading  parameter. 

CC4  boundary  conditions  were  assumed  for  the  buckling  solution  (re¬ 
call  that  prebuckling  displacements  were  neglected),  which  means 
u°(O,0)  =  v° (0,0)  =  w(0,9)  =  w,x(O,0) 

(A7) 

-  u°(L,0)  =  v°(L,0)  =  w(L,8)  =  w,x(L,0)  =0 
where  L  denotes  cylinder  length.  A  solution  to  Eqs.  (Al)  -  (A3)  is  a 
Fourier  series. 

The  displacementsw  becomes 

w  ■  sin  n8  T  F  sin 

**  mn  I . 


,  A  r  „  miTX 

+  cos  n0  )  G  cos  — — 

mn  L 

m*0 
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where  n  and  m  are  1/2  sine  or  cosine  waves  in  the  circumferential  and 
axial  directions  respectively.  The  displacement  is  continuous  in  the  6 
direction  so  there  is  no  need  to  sum  over  the  terms  containing  9 .  The 
displacements  u°  and  v°  are  of  similar  form  as  w. 

Taking  the  equilibrium  equations  (Al)  -  (A3),  along  with  the  dis¬ 
placement  functions  for  u°,  v°  and  w°,  and  the  boundary  conditions,  Eq. 
(A7),  yields  in  general,  two  6x6  sets  of  equations  for  each  value  of  n. 
Setting  up  a  coefficient  matrix  for  the  6x6  system  and  setting  it 
equal  to  zero  will  yield  values  of  n.  The  buckling  load  corresponds  Lo 

the  value  of  n  which  yields  the  lowest  value  of  N  . 

o 

Since  the  Fourier  procedure  can  become  very  cumbersome,  a  method  to 
determine  a  quick,  estimate  is  often  used.  Whitney  and  Sun  used  this 
approximation  to  compare  to  the  Fourier  analysis.  The  approximation  as¬ 
sumes  for  the  solution  to  Eqs.  (Al)  -  (A3)  the  following  form 
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Appendix  B 

User  Problems  with  STAGS 

In  operating  STAGS  during  this  thesis,  a  number  of  problems 
occurred  that  were  caused  by  certain  input  parameters  and  models.  As 
a  result  of  the  problems,  many  weeks  were  taken  from  research  of 
composite  shells  to  researching  the  operation  of  STAGS.  The  cases 
that  are  described  in  this  appendix  should  aid  any  follow-up  work  to 
this  project. 

Shift  in  Eigenvalue 

An  initial  eigenvalue  shift  (shift  on  F-3  card)  of  1.  worked  for 
all  cases  in  this  thesis.  Values  of  0.  and  3.  were  both  unsuccessful 
as  they  had  the  eigenvalue  converge  to  a  negative  root.  Also,  E1GA= 
EIGB=NEIG=1,  was  also  successful. 

Rigid  Body  Displacements 

Loading  of  the  full  shell  with  torsion  at  both  the  top  edge  and 
bottom  edge  (in  opposite  directions!)  should  not  generate  displacements 
of  a  rigid  body.  Unfortunately,  it  does.  The  way  around  this  situ¬ 
ation  is  to  hold  at  least  one  node  on  an  edge  fixed  in  the  circum¬ 
ferential  direction.  It  was  discovered  that  holding  all  the  nodes  of 
one  edge  fixed  in  the  circumferential  direction  made  for  an  "easier  to 
follow"  displacement  vector.  It  is  also  easier  to  input.  Simply 
apply  the  shear  to  one  edge  with  an  L-2  card  (STAGS-C)  and  let  V=0 
on  the  K-2  and  K-4  card  for  the  other  edge.  This  will  change  the 
boundary  condition  on  that  edge  but  not  the  edge  with  the  applied 
shear  and  it  should  run  smoothly. 
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Ill-Conditioned  Variables 

The  ill-conditioned  variables  are  due  to  numerical  instability. 

The  solution  is  to  simply  change  your  mesh  size.  Increasing  the  mesh 
size  by  around  15%  corrected  the  problem  for  the  case  mentioned  in  this 
thesis. 

Additional  Information 

It  is  important  to  remember  that  when  applying  a  load  to  two  edges 
of  a  full  shell  (axial  compression,  torsion,  etc.),  an  additional  L-2 
card  is  required.  It  is  not  an  additional  load  case  if  the  loads  on 
the  two  edges  are  not  independent  so  N  LOAD  =  1  on  the  1-1  card. 
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TABLE  1 


AXIAL  COMPRESSION  -  =  7.0 
K 


CCl 


CC2 


N„/N„ 


N 

X 

v  xss 

N 

X 

A  SS 

0° 

56933 

.7591 

63780 

.8504 

15° 

82311 

.9407 

101509 

1.1601 

30° 

102247 

.8891 

116679 

1.0146 

45° 

128313 

.9684 

128379 

.9689 

60° 

131515 

.9070 

131515 

.9070 

90° 

66293 

.8839 

66293 

.8839 

0 


CC3 


CC4 


N 


x 


N  /N 
x  x 


SS 


N 

x 


N  /N 
x  x 


SS 


0 

56978 

15° 

82320 

30° 

102787 

45° 

128340 

60° 

131515 

90° 

66293 

.7597 

.9408 

.8938 

.9686 

.9070 

.8893 


63810 

102944 

116725 

128379 

131515 

66293 


.8508 

1.1765 

1.0150 

.9689 

.9070 

.8839 
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TABLE  2 


AXIAL  COMPRESSION  -  =  12 
R 


0 

CC1 

CC2 

Nx 

N  /N 

N  /N 

x  XSS 

Nx 

X  X„„ 

SS 

0° 

58561 

.7808 

58733 

.7831 

15° 

77901 

.8903 

95043 

1.0862 

30° 

113/00 

.9887 

113724 

.9889 

45° 

„  o 

127306 

.9608 

127306 

.9608 

60 

131385 

.9061 

131385 

.9061 

90° 

66173 

.8823 

66173 

.8823 

0 

CC3 

CC4 

Nx 

N  /N 
x  XSS 

N 

X 

N  /N 
x  x„„ 
SS 

0° 

59025 

.7870 

59040 

.7872 

15° 

80448 

.9194 

98105 

1.1212 

30° 

113712 

.9888 

113735 

.9890 

45° 

„  o 

127319 

.9609 

127319 

.9609 

60° 

o 

131385 

.9061 

131385 

.9061 

90 

66173 

.8823 

66173 

.8823 
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TABLE  3 


TORSION  L/R  =  7 


CC1  CC2 


Nx0 

N  /N 
x0  x6ss 

Nx0 

N  /N 

X0  X0 

15° 

19533 

1.1490 

20349 

1.1970 

45° 

46308 

1.0896 

48043 

1.1305 

60° 

58939 

1.1671 

60591 

1.1998 

75° 

56916 

1.0841 

60434 

1.1511 

90° 

49866 

1.1206 

50305 

1.1305 

CC3  CC4 


N  /N 

N  /N  n 

N 

X0  X0 

N  Q 

X0  X0C 

X0 

SS 

X0 

8° 

16025 

1.0865 

16375 

1.1102 

15° 

19684 

1.1579 

20489 

1.2053 

45° 

46650 

1.0977 

48340 

1.1374 

60° 

59997 

1.1881 

60757 

1.2031 

75° 

57474 

1.0947 

60967 

1.1613 

90° 

49962 

1.1227 

50376 

1.1320 
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TABLE  4 


TORSION 

L/R  =  12 

0 

CC3 

CC4 

NX0 

Nx0/Nxess 

Nx0 

N  «/N  . 
x0  x0( 

8° 

13070 

.8861 

13147 

.8913 

15° 

15312 

.9007 

15912 

.9360 

45° 

32950 

.7753 

37489 

.8821 

60° 

42344 

.8385 

45415 

.8993 

75° 

47859 

.9116 

48799 

.9295 

90° 

47735 

1.0727 

48189 

1.0829 
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